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The time-dependent density functional based tight-binding (TD-DFTB) approach is generalized 
to account for fractional occupations. In addition, an on-site correction leads to marked qualitative 
and quantitative improvements over the original method. Especially, the known failure of TD-DFTB 
for the description of a — > it* and n — > tv* excitations is overcome. Benchmark calculations on a 
large set of organic molecules also indicate a better description of triplet states. The accuracy of the 
revised TD-DFTB method is found to be similar to first principles TD-DFT calculations at a highly 
reduced computational cost. As a side issue, we also discuss the generalization of the TD-DFTB 
method to spin-polarized systems. In contrast to an earlier study [Trani et al., JCTC 7 3304 (2011)], 
we obtain a formalism that is fully consistent with the use of local exchange-correlation functionals 
in the ground state DFTB method. 
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I. INTRODUCTION 

During the last years, Density Functional Theory (DFT) has been one of the most utilized tools for the description 
of ground-state properties of a wide variety of molecular systems that range from small molecules to large periodic 
materials. While it lacks the accuracy typical of correlated wavefunction-based methods, it goes beyond Hartree-Fock 
(HF) as electron correlation is incorporated in a self-consistent-field (SCF) fashion. Thus, DFT has turned out to 
be a good compromise between accuracy and computational cost; affordable to study hundreds-of-atoms systems on 
most of the current work stations with fairly good precision. The field of application of this method was subsequently 
extended to the study of excited states properties with the development of time dependent density functional theory 
(TD-DFT) P This method has become the de facto standard for the computation of optical properties for molecules 
with several tens of atoms. Also the limitations of TD-DFT are now well documented in the literature (see e.g!^), 
which allows researchers to judge a priori whether a certain class of exchange-correlation functionals is sufficient for 
the predictive simulation of the problem at hand. 

Still, many applications in photochemistry and nanophysics are not easily tractable in the current methodological 
frame. As an example, quantum molecular dynamics in the excited state require the evaluation of energies and forces 
at a large number of points along the trajectory. Also, the investigation of extended nanostructures like nanowires 
and -tubes with intrinsic defects or surface modifications can not be reliably performed with small models. For these 
kind of problems an approximate TD-DFT method might be advantageous. Such a scheme is the time dependent 
density functional based tight binding method (TD-DFTB) ™ In the TD-DFTB method additional approximations 
beyond the choice of a given exchange-correlation functional are performed to enhance the numerical efficiency. These 
are mostly the neglect and simplification of two-electron integrals occurring in the linear response formulation of 
TD-DFT. In contrast to HF-based semi-empirical methods (like INDO/S or CIS-PM3), TD-DFTB approximates a 
reference method that already covers electronic correlation and does not include free parameters that are adjusted to 
reproduce experimental data. 

After the original linear response implementation^ TD-DFTB was extended in a number of different directions. 
Analytical excited state gradients have been derived^ as well as a real time propagation of Kohn-Sham (KS) orbitals 
that enables the computation of optical spectra using order-N algor ithm s!^ Several groups worked on non-adiabatic 
molecular dynamics simulations in the EhrenfeslP or surface hopping^l approach. More recent developments include 
a TD-DFTB approach for open boundary conditions in the field of quantum transporters well as an implementation 
for open shell systems.^ A detailed review that summarizes advantages and limits of the method is also available.^ 

One of these limits is the rather poor description of a — > n* and n — > tt* excitations in many chromophores. In 
TD-DFTB these transitions show vanishing oscillator strengths although higher levels of theory clearly predict a finite 
(albeit small) transition probability. Moreover, singlet-triplet gaps for these kind of excitations are exactly zero in TD- 
DFTB, which is again in disagreement with more accurate methods. Although n — > tt* transitions usually dominate 
the absorption spectrum, such low-lying excitations may play a significant role in the luminescence properties that 
are of key importance in many technological applications. 

In this paper, we present extensions of the TD-DFTB method in order to improve its description of the absorption 
spectra of molecules. In the first section, the spin-unrestricted DFTB formulation for the ground state is briefly 
reviewed. Afterwards, we generalize the TD-DFTB approach for the study of open-shell systems with fractional 
occupancy. In |IV| a refinement to overcome important deficiencies within the method is formulated. In order to 
highlight the qualitative improvement due to these extensions, we report results for selected diatomic molecules in 
|V} Additionally, a comparison between results obtained with the proposed formalism and the original TD-DFTB 
approach for a large set of benchmark molecules is presented. Our findings are further compared to TD-DFT, the 
best theoretical estimates from the literature and experimental observations. 

II. SPIN-UNRESTRICTED DFTB FOR THE GROUND STATE 

This section contains a brief summary of the spin-unrestricted DFTB method in order to establish the required 
notation for the later parts. More details on the derivation and application of the method are given in the original 
article by Kohler and co-workers.^ 

DFTB is based on a expansion of the DFT-KS energy up to the second order in the charge density fluctuations, 
Sp = Sp-f + 5pi, around a reference density, po = [pof,Poi\- The latter is given by a superposition of pre-calculated 
densities for neutral spin-unpolarized atoms. The DFTB energy reads: 
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+ 5 P^ v ) f^Pr [Po](r, r') ^ r (r') drdr' + £ rcp , (1) 

where the rii CT denote occupation numbers, E^-cp is a sum of pair potentials independent of the electronic density,^ 
and / hxc = f h + ,r c with 



/ h (r,r') = - - - 
|r — r' 

S 2 E xc [p] 

5p a {v)5p T (v') ! 



KprWM = Sn M ?r, W v ( 2 ) 



denote the Coulomb and exchange-correlation kernel, respectively. 

In the first (zeroth-order) term of [l] H° is the KS Hamiltonian evaluated at po and the sum runs over the spin 
variables a =f, i and the KS orbitals ipia- For the succeeding formulation, the Roman indices s, t, u, v, denote KS 
orbitals, whereas Greek indices p, v, k, X will denote atomic orbitals (AO). Let us also abbreviate a general two-point 
integral over a kernel g(v, r') in the following form: 



(f\g\h) = JI f(v)g{v, v')h(r') dvdv'. (3) 

For atomic orbital products /(r) = M (r)</)„(r) and h(v') = cj) K {v')(f)x{v') the shorthand {pv\g\n\) will be also used. 
Expanding the KS orbitals (which we choose to be real) into a suitable set of such localized atom-centered AO, 
V'io- = c JIi0M' t ne second term in |TJ labeled in the following, can be written as 

E(2) ^EE AP ^ (H/£l \Po]\**) &Kx, (4) 

(TT flUKX 

where AP^„ = PZ V — PJJjJ 7 denotes the AO density matrix of the difference density 8p a (r) , with 



P° 



EZ c p s 



per a p 
1 st c vti r i 



■O.cr 



'^ILS ± St '-V 



(5) 



Here P° t = (ty\a t<7 a ser \ i f>) are the MO density matrix elements, ^ being the ground state KS determinant. The term 
P s f designates the MO density matrix of the reference system. 

To evaluate the appearing integrals, the Mulliken approximation is applied. This amounts to set fy^v ~ 
|<S/ii/(|<^Vl 2 + l^l 2 )' using the known AO overlap integrals S^. Thus, the general four-center integrals are writ- 
ten in terms of two-center integrals, 

H/^ c t ni«a)«^ M!/ 5 kA ]t Yl M/p h ;p>o]i/^) • ( 6 ) 

By substituting [6] into [4j the second order contribution to the energy reads, 



Ei2) = ^EE A ^.W^WH A ^. ( 7 ) 

(TT 11V 

where P° are the elements of the dual density matrix defined as 



P = - (P ■ S + S ■ P) . 



(8) 
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By using [5] the dual density matrix can be expressed with respect to the MO density matrix as 

K = J2 P ^ P sf (9) 

St 

The matrix Pj^J 7 represents the Kohn-Sham transition density for an excitation from orbital s to t, and is defined as 
follows: 



Pflt ~ ^ { c tis^i/t + c )it^vs + c vs^iit + c vt^fj,s) ! c s — c s ■ S. (10) 

A further simplification to [7] is obtained by spherical averaging over AO products. This ensures that the final total 
energy expression is invariant with respect to arbitrary rotations of the molecular frame. To this end, the functions 
Fai are introduced 

m—l 

F Al (\r-R A \) = WTl £ |^(r-R x )| a , (11) 

m=—l 

where I and m denote the angular momentum and magnetic quantum number of AO (f>a, centered on atom A. We 
then have for p, = {Aim}, v — {Bl'm 1 }: 

(M/p;, C >o]M « {F M \f^ p M\ F Bv) = ?M,Bl> (12) 

introducing the shorthand notation T. 

The traditional DFTB energy expression is now obtained after transformation from the set {pf,p±} to the total 
density p = p^ + Pi and magnetization m = p^ — p±. By this change of variables, the exchange-correlation kernel can 
be split and one arrives at 



^m,bv = 7ai,bi> + S<t5 t SabWai,i>, (13) 

where S a = 2tf CTt - 1 and the parameters Jai,BI< = {F A i\f^ c [p ]\F B v) and Wai,v = {F A i\f™ m [p ]\F AV ). Please 
note that the constants W A i^i depend only on the exchange-correlation kernel, but not on the long-range Coulomb 
interaction. Moreover, as the reference density po is built from neutral spin-unpolarized atomic densities, there are 
no integrals which involve mixed derivatives of the exchange-correlation energy with respect to both density and 
magnetization. The parameters Jju.bi' and W A iy are known in the DFTB method as the 7-functional and spin 
coupling constants, respectivelyP^l The spin coupling constants are treated as strictly on-site parameters, whereas 7 
is calculated for every atom pair using an interpolation formula, that depends on the distance between the atoms 
A and B and the atomic Hubbard-like parameters r y A i. A i and jbi,bi- Traditionally, the latter are not computed 



directly from the integral 12 but from total energy derivatives according to 7a;, ai = S 2 E/8n 2 . In this expression, E 
denotes the DFT total energy of atom A and n refers to the occupation of the shell with angular momentum I. The 
derivative is then evaluated numerically by full self-consistent field calculations at perturbed occupations. Due to 
orbital relaxation, parameters obtained in this way are roughly 10-20 % smaller than the ones from a direct integral 
evaluation. This point will become important later. 

Finally, by substituting [9] in [7] while using [12] the second order energy term can be expressed as 

E(2) ^P>mTai,bv<Iw^PL, (14) 

err stuv ABU' 

where q A i = X^eA / are the angular-momentum-resolved transition charges for atom A. For later reference, we 
next define a matrix K: 

J^sta,uvT — "Al 1 Al,BVWl' ■ \ L0 ) 

ABU' 

Using this abbreviation, the DFTB total energy can then be written in the following simple form, 
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E = Y.H H °staP!t + \Y.ll^ P stKsta^P T uv +E rcv . (16) 

a st err stuv 

By applying the variational principle to the energy functional, the DFTB KS equations are obtained 

H sta - e S(T S st = 0, Vs,t,a (17) 

where the KS Hamiltonian is expressed as 

□ El 

H «° = gp^ = H° sta + J2Yl K st ^ UVT AP: v . (18) 
s ^ r uv 

These equations are subsequently transformed into a set of algebraic equations by expanding the KS orbitals into 
the AO basis, and are solved self-consistently. For the converged ground state density we have P° t = n sa 8 s t- Using 
also the identity J2st c 1isPsi' T c Zt = n %o^v (V//, ^, a - ), n® being the occupation numbers for the reference atoms, it is 
straightforward to recover the expression for the spin-unrestricted DFTB Hamiltonian given ir l 12 * 13 [ 

III. SPIN-UNRESTRICTED TD-DFTB 

Once the KS orbitals ipi a and energies Ci a are obtained, a linear density-response treatment directly applies as a 
natural extension to DFTB. Following Casidap31the excitation energies ujj and eigenvectors Fj can be calculated by 
solving the eigenvalue problem 

OF, = ujFj, (19) 

where the response matrix CI is defined as 

fi = S- 1 / 2 (A + B)S- 1 / 2 (20) 

S = -C(A-B) _1 C. (21) 
The matrices A, B and C are defined according to 





SijSabSarUJjbr 




Tlj T T^br 




^iaa^bjr 




_ SijSabSfjr 


71 jf Tlbr 



K. 



iaa,jbT 



(22) 

where u^jbr = e br ~ e jr arL d the indexes a, b stand for KS orbitals such that n i(T > n aa and rij T > n^. We explicitly 
allow for fractional occupations at this point in order to be able to simulate also metallic or near-metallic systems at 
elevated electronic temperature. 

The coupling matrix elements Ki aa ^i, T are generally defined as the derivative of the KS Hamiltonian with res pect 



to the density matrix elements. For the special case of DFTB (see 18 1, this leads to the matrix K defined in 
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This quantity depends only on the 7 and W parameters, as well as on the Mulliken transitions charges obtained from 
the previous ground state DFTB calculation. 



G 



An important feature of the coupling matrix for local or semi-local exchange-correlation functional is its invariance 
with respect to the permutation of any connected (real orbital) indices (e.g. Ki aa ji, T = Ki aai bj T — K a i a- jb T )- This 
symmetry does not hold for functionals involving Hartree-Fock exchange.^ In contrast to an earlier derivation,^ W e 
find that the DFTB coupling matrix is in fact symmetric. This can be traced back to the transition density matrix, 



10 that evidently has this property. Since the ground state DFTB method is derived as an approximation to local 
or semi-local DFT only2^, this is an expected and necessary property and implies that the orbital rotation Hessian 
A — B becomes strictly diagonal. Thus, the expression for the response matrix elements is conveniently simplified to 
read 

Qiaajbr = SijSabSar^'jbT + 2-\/ (n,i a — n aa ) WiaaKiaaJbr "J (%t — n br) ^jbr- (24) 

It is worth formulating the method for closed shell systems as a particular case, in order to make contact with the 
original derivation of TD-DFTB . In this special case the Mulliken transition charges have the property q % %} = q™i = 
q^i- If, in addition, the dependence of the 7- functional and W constants on the angular momentum is neglected, the 
coupling matrix simplifies to 

Kiaajbr = q A + 5 a 5 T 5 AB W A ) q$, (25) 

AB 

with q™ = q At , in full agreement with the expression derived previously for spin-unpolarized densitiesP 

Once the eigenvalue problem is solved within TD-DFTB, the oscillator strength related to excitation / can be 
calculated as 



9 3 

/<= E 

k = l 



(26) 



where denotes the fc-th component of the position operator. The transition-dipole matrix elements are also subjected 
to a Mulliken approximation (Ra is the position of atom A): 

(^ CT |r|^ OCT )«^RA<dr. (27) 



IV. BEYOND THE MULLIKEN APPROXIMATION 

While in general the Mulliken approximation accounts, at least approximately, for the differential overlap of atomic 
orbitals, there is an important exception. If orbitals M and 0„ with /1 / y reside on the same atom, their product 
vanishes since the overlap matrix reduces to the identity. As a consequence, the transition charges are underestimated 
or often vanish identically for excitations that feature localized MO, e.g. n — > n* transitions. This in turn means that 
the coupling matrix vanishes and hence no correction of the ground state KS energy differences occurs in the linear 
response treatment. More importantly, the singlet-triplet gap is zero for such excitations. 

A next level of approximation demands therefore the evaluation of all one-center integrals of the exchange type, 
i.e., (fJ-v\f^Z[po]\fJ>v) with (i / 1/. This is how Pople et al. proceeded in the development of the intermediate 
neglect of differential overlap model (INDO)P^ to overcome the deficiencies encountered within the complete neglect 
of differential overlap method (CNDO)P^ Under the INDO model, the one-center, two-electron integrals are fitted to 
atomic spectroscopic data. In this work, the corresponding one-center integrals are calculated by numerical integration. 

After inclusion of every onsite exchange-like integral within the DFTB formalism (see [A]) , the second-order energy 
takes the new form: 
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E(2) = J EE E Ai* (H/^boiMAi* 



2 



EE E A ^ H^WH ap;, 

A^B 

o E E E E a^ H^WH ap;„ (28) 
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where the onsite integrals appear in the second term and the previous energy expression |7| has been split into on-site 
(first term) and off-sit e (th ird term) c ont ributions. 

By substituting [9] in [28] while using 12 for the off-site component, E^> can be finally written as 



E(2) = \ E E &p: t K s t„,u V rAP: v , (29) 



2 

<7T StUV 

with the refined coupling matrix K being expressed as 

w -EE ^Hj^whc 

A [iveA 

+ 2 E E p *H/jr>]iH^ T 

A fiveA 
A^B 

+ E <fM?M,BV<lm- (30) 

ABIV 

With inclusion of the exchange integrals, the spherical averaging over AO products described in |12| will in general 
not lead to an expression which is invariant under a rotation of the atomic axes. This has been discussed by Figeys 
et. alpS who analyzed the rotational invariance (RI) of the INDO model. Taking p-type orbitals as an example, the 
following identity must hold to preserve RI: 

(j>p\f?: c PT [po]\p P ) - (H/ P ^p>oW) - 2(pp'\f£ c pT [p ]\pp'), Vp,p' . (31) 

In the original DFTB formulation this requirement is fulfilled, as both integrals on the left-hand side of [31] are 
approximated to have the same value, T'jf A , while the integral on the right-hand side is neglected. To retain RI 



within the present scheme, one could evaluate all on-site integrals exactly so that 31 holds automatically. As detailed in 
ITU the corresponding averaged T-parameters would then differ from the ones usually employed in ground state DFTB 
simulations, which are obtained from total energy derivatives. Instead, we use [Ml to calculate all required parameters, 
taking the exact exchange integrals and the traditional T-parameters as inputP^ In this way the modifications of the 
original method are kept as small as possible, while RI is still exactly preserved. 

The eqs [29] and [30] represent a correction to the conventional DFTB energy expression, whose implication for various 
ground state properties certainly warrants a deeper investigation. Since we are mainly interested in an improvement 
of the TD-DFTB scheme at this point, we keep with the traditional DFTB scheme for the ground state and employ 



the modified coupling matrix only in the response part of the calculation, i.e. in 22 

In a similar manner as the refinement of the coupling matrix was achieved, the approximation for the dipole matrix 
elements eq |27[ and hence the oscillator strength, can be improved by including all non- vanishing one-center dipole 
integrals (see |B[), 

(vwifivw) = J2 K ^A a + E E p ;r<Mi%). (32) 

A A /use A 

According to the dipole selection rules only s-p and p-d dipole integrals are non-zero. Among the s-p integrals, only 
those of the type D sp — (s\rk\pk) do not vanish, all of them being equal. With regards to the p-d integrals, eleven of 
them are non-vanishing: 
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(p y \x\d xy ) = (p x \y\d xy ) = (p y \z\d yz ) = (jp z \y\d y z ) = (p z \x\d xz ) = (p x | z I d xz ) 

(P X \x\d X 2_y2) = -(Py\y\d X 2_y2) 

(p y \y\d z 2) = (p x \x\d z 2) 

( Pz \z\d z 2). (33) 

We would like to stress that the required additional atomic integrals in our correction to the coupling matrix and 
oscillator strength are neither freely adjustable parameters nor fitted to experiment. Instead, they are calculated 
numerically with a special-purpose DFT code that was recently developed in our group. 

Like in earlier studies with the DFTB method, the Perdew-Burke-Ernzerhof (PBE}E0 functional was used in the 
calculation of Hamiltonian and overlap matrices, as well as the computation of reference densities and basis functions.^ 
This also holds for the new parameters, which are calculated once for every atom type, stored in a file, and read when 
needed during the calculation. 

The presented extensions of the TD-DFTB method have been implemented in a development version of the DFTB+ 
codeP 



Dpd — 
Dpd — 
D'pd = 



V. RESULTS 



A. Diatomic molecules 



To show the performance of our corrections to the coupling and dipole matrices (hereafter referred to as the onsite 
correction) , we calculated the low- lying vertical excitation energies and their corresponding oscillator strengths of three 
diatomic molecules for which the improvements over traditional TD-DFTB are especially noticeable. It is important 
to point out that only valence excited states can be treated within TD-DFTB due to the employed minimal basis set. 
[T] shows the excitation energies of NO, N2 and O2 calculated within both the corrected and the original formulation of 
TD-DFTB. For NO and O2, spin-unrestricted TD-DFTB calculations have been performed, where doublet and triplet 
ground states have been considered, respectively. In order to identify the excited state multiplicity of the open-shell 
systems, the expectation value of the square of the total spin operator, (S* 2 ), is evaluated. We apply the expression 
derived in RefP3 for unrestricted configuration interaction singles (UCIS) wave functions for this purpose. In |TJ a 
multiplicity is only assigned to those excited states with low spin contamination. This covers the most important 
excitations in the absorption spectrum. 

As a reference, we computed vertical excitation energies to valence states of these compounds by using TD-DFT as 
implemented in the TURBOMOLE package^. The PBE exchange-correlation functional as well as triple-zeta plus 
polarization (TZP) basis set has been used. All ground state geometries were previously optimized at the cor respo nding 
level of theory. Some experimental findings taken from the literature are also included for compariso n 25 ^ 26 ! . The 
oscillator strength for each excitation is indicated to the right of the corresponding excitation energy and in the first 
column of the table the symmetry and type of the transition is specified. 

These molecules have as a common feature that they all possess low-lying a — > 7r* excitations playing an important 
role in their absorption spectra. As stated above, a wrong description of these transitions is a known issue in current 
TD-DFTB. Below, we identify yet another failure in the description of some ir — > it* transitions of these compounds. 

NO belongs to the symmetry point group Coot, for which LI and E + transitions are electric dipole allowed. However, 
TD-DFTB describes the former as forbidden. This is due to the mentioned vanishing of the corresponding transition 
charge which leads to an equality of the KS energy difference wks (termed ujjbr m |22|) and the excited state energy 
uji. Within the refined formulation this failure is successfully overcome as shown in III This improvement is specially 
important for the second and fourth LI transitions with oscillator strengths of 0.01 and 0.24, respectively, which are in 
agreement with the TD-DFT values of 0.02 and 0.38. Our correction is in this case, essential for providing a correct 
absorption spectrum (see[TJ where traditional TD-DFTB is able to describe only the S + peak. 

The onsite correction also improves the description of II transitions quantitatively For example, the first II excita- 
tion energy is clearly overestimated with respect to first-principle results. When using our correction the overestimation 
is reduced by almost 0.5 eV. Oppositely, the second II excitation energy is strongly underestimated within traditional 
TD-DFTB whereas the corrected calculations return a value (8.33 eV) close to that from TD-DFT (8.61 eV). The 
S + (ir —> it*) transitions are also found to be better described within our correction, although the first (second) S + 
excitation energy is still significantly overestimated (underestimated). The £~ (ir — > tt*) excitations energies are on 
the other hand unaffected by our correction. For these transitions no shift from their KS energy differences are seen, 
which agrees with first-principle observations. More importantly, the onsite correction rectifies the wrong degeneracy 
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Molecule/Trans. TD-DFT TD-DFTB Exp. 
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0.00 


7.91 


6.80 


0.00 


6.79 


0.00 


6.79 




3 V + 


- 




9.00 


0.18 


6.90 


8.36 


0.32 


8.21 


0.24 


6.35 


~8.6 


3 n 

1 ' a 






14.85 


0.18 


14.16 


15.35 


0.18 


14.52 


0.00 


14.52 





TABLE I: Comparison of vertical excitation energies Wi and oscillator strengths fi for TD-DFT with TZP basis set, the 
traditional TD-DFTB method (old) and TD-DFTB with onsite correction (new). The PBE functional is used throughout. 
o;ks denotes the KS orbital energy difference corresponding to the most dominant single particle transition in the many body 
wavefunction, as discussed by Casida in RefP^. Experimental data for N 2 and 2 were taken from RefP^ and inferred from 
the potential energy curves in RefPP, respectively. Oscillator strengths are only provided for excitations that are not trivially 
spin- forbidden. All energies are expressed in eV. 

of the transitions £~ and A, which we identify as another important qualitative issue within traditional TD-DFTB. 
The correction also reduces the spin contamination of the doublet 2 II states. 

The electric-dipole-allowed transitions 1 H U and 3 IT U of the homonuclear molecules N2 and O2 (point group 
Dooft)) respectively, are neither correctly described by traditional TD-DFTB as shown in [i] When applying the onsite 
correction, these excitations become allowed with oscillator strengths in agreement with ab-initio results. However, it 
is necessary to indicate that the excitation energy of the latter transition is somewhat overestimated, being in better 
agreement within the non-corrected formalism. In contrast, the onsite correction greatly improves the agreement of 
the excitation energy of N2 with the TD-DFT result. This transition and the 3 H U state are degenerate according 
to traditional results whereas such degeneracy is broken at the corrected TD-DFTB level. This degeneracy breaking 
is in total correspondence with the results obtained at a higher level of theory. 

The original formalism also predicts the degeneracy of the singlet and triplet £~ and A u states of N2 with excitation 
energy equal to 9.01 eV. This is however only partially confirmed by the ab-initio results and in total disagreement 
with the experimental findings. According to TD-DFT, the triplet and singlet E~ states are degenerate with corre- 
sponding excitation energy of 9.60 eV but the triplet and singlet A u degeneracy is not observed. On the other hand, 



TD-DFT O TD-DFTB (new) O TD-DFTB (old) 



10 




FIG. 1: Absorption spectrum of nitric oxide as obtained with full TD-DFT (PBE/TZP), traditional TD-DFTB (old) and 
TD-DFTB with on-site corrections (new). 



experiments report nearly degenerate S~ states with excitation energies of 9.67 and 9.92 eV for the triplet and singlet 
states, re specti vely, in contrast with TD-DFT findings. This apparent failure of TD-DFT has been also reproduced 
elsewhere! 27 * 28 ^. In a recent letter it was shown that excited states like N 2 E~ cannot be described by linear-response 
TD-DFT as their corresponding excitation energies do not correspond to poles of the response functiorP3. TD-DFTB 
as an approximation to TD-DFT unavoidably inherits this issue and our correction is unable to fix it. However, 
it does break the wrong degeneracy of the triplet and singlet A u states. For O2, a similar degeneracy breaking of 
the transitions S~ and A u within the corrected TD-DFTB can be observed. This is again in total agreement with 
first-principle results. 



B. General benchmark 



To assess the general performance of the corrected TD-DFTB method in terms of vertical excitation energies we 
have chosen a large benchmark set defined elsewhere^ and recently used for the validation of TD-DFTB^4 This 
set covers n — > tt* , n — > n* and a — > 7r* excitations for 28 organic compounds comprising unsaturated aliphatic 
hydrocarbons (group 1), aromatic hydrocarbons and heterocycles (group 2), carbonyl compounds (group 3) and 
nucleobases (group 4), and intends to embrace the most important chromophores in organic photochemistry. For 
comparison we have calculated singlet and triplet vertical excitation energies at the TD-DFT (PBE, TZP) level. To 
disentangle the accuracy of the TD-DFTB approximation itself from the quality of DFTB ground state geometries, 
the same optimized structures were used along the benchmark. These geometries were previously obtained at a MP2 
levePSJ In the Supporting Information, an extended table containing all our calculation results can be consulted. For 
further comparison, we also provide the theoretical best estimates (TBE) for this benchmark set, calculated at the 
CASPT2 level and reported by Schreiber and co-workersP^l. Reported experimental values for some vertical excitation 
energies are additionally given. We also include the KS energy difference of the corresponding dominant single-particle 
transition, which is useful for the analysis of the relative displacement of the excitation energies with respect to the 
KS energy difference as compared to ab-initio TD-DFT. 

In[lI|and|III[ we present a statistical analysis of the collected data for the groups of compounds previously mentioned, 



for singlet-to-singlet and singlet-to-triplet transitions, respectively. In IV the global statistics for the benchmark set 
is summarized. Here the analysis has been also split into triplet and singlet excitations. The mean signed difference 
(MSD) as well as the root-mean-square deviation (RMS) of the TD-DFTB excitation energies with respect to the 
experiment, the TBE and TD-DFT are reported. Statistics on ab-initio TD-DFT results are also provided to indicate 
its degree of correspondence with higher level of theory and experiment. 

To assess the validity of our approximation, we focus on the comparison with the TD-DFT (PBE) data set. It is 
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count TD-DFT TD-DFTB 









new 


old 


Group 1 










MSD (TD-DFT) 


13 




0.15 


0.07 


(TBE) 


13 


-0.52 


-0.37 


-0.45 


(Exp.) 


12 


-0.35 


-0.20 


-0.26 


RMS (TD-DFT) 


13 




0.28 


0.25 


(TBE) 


13 


0.62 


0.47 


0.53 


(Exp.) 


12 


0.56 


0.39 


0.46 


Group 2 










MSD (TD-DFT) 


53 




0.23 


0.12 


(TBE) 


53 


-0.35 


-0.12 


-0.23 


(Exp.) 


12 


-0.06 


0.13 


0.02 


RMS (TD-DFT) 


53 




0.41 


0.36 


(TBE) 


53 


0.54 


0.37 


0.43 


(Exp.) 


42 


0.39 


0.35 


0.34 


Group 3 










MSD (TD-DFT) 


19 




0.56 


0.37 


(TBE) 


19 


-0.81 


-0.25 


-0.44 


(Exp.) 


13 


-0.67 


0.07 


-0.08 


RMS (TD-DFT) 


19 




1.06 


0.93 


(TBE) 


19 


0.96 


0.88 


0.88 


(Exp.) 


13 


0.88 


0.70 


0.64 


Group 4 










MSD (TD-DFT) 


19 




0.09 


0.03 


(TBE) 


19 


-0.84 


-0.75 


-0.81 


(Exp.) 


13 


-0.51 


-0.37 


-0.29 


RMS (TD-DFT) 


19 




0.32 


0.30 


(TBE) 


19 


0.89 


0.91 


0.96 


(Exp.) 


13 


0.63 


0.61 


0.64 



TABLE II: Mean signed difference (MSD) and root-mean-square deviation (RMS) of singlet-singlet vertical excitation energies 
(in eV) calculated within TD-DFTB and TD-DFT (along the rows) with respect to TD-DFT, TBE and experiment (along the 
columns) where possible. 

important to recall that TD-DFTB parameters are calculated at this level of theory and the main aim of our approach 
is to improve its agreement with respect to the TD-DFT description. Since TD-DFT at the PBE level is of course 
an approximation itself, we are also interested in examining the accuracy of our method compared to a higher level 
of theory and experiment. It is worth mentioning that the subset of singlet-triplet excitations for which experimental 
observations are provided is around 50% smaller than the original set, but we consider it still suitable to perform a 
statistical analysis. The TBEs are available for the complete benchmark set on the other hand and they are fairly 
close to their corresponding experimental values, with a RMS error of 0.24 eV for triplet states and 0.38 eV for 
singlet states. On average, both singlet and triplet TBEs are slightly overestimated (MSD = 0.15 e V and 0.20 eV, 
respectively) which should be taken into account in the further analysis. Consistent with earlier studies EU3H TD-DFT 
excitation energies at the PBE level appear to be somewhat underestimated with regard to experimental findings as 
a general trend (MSD = -0.28 eV and -0.27 eV for singlet and triplet states, respectively). 

It can bee seen from |IV| that with respect to TD-DFT, traditional TD-DFTB performs better for singlet-singlet 
than for singlet-triplet excitations. It should be pointed out however, that triplet states are in better agreement with 
the TBE values as it has been previously noticed 1 ^. Nevertheless, compared to the collected experimental findings, 
singlet states also appear to be better described. In fact, |IV| shows a significant overestimation of triplet state energies 
compared to experiment (MSD = 0.50 eV) within traditional TD-DFTB. One of the main effects of our correction 
is the significant improvement of singlet-triplet excitation energies taking TD-DFT and experimental values as a 
reference. Within the refined formulation, the RMS error is reduced by 0.3 eV with respect to both TD-DFT and 
experiment. More importantly, the RMS error with respect to the latter (0.33 eV) is slightly lower than that for 
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count TD-DFT TD-DFTB 









new 


old 


Group 1 










MSD (TD-DFT) 


13 




0.21 


0.71 


(TBE) 


13 


-0.32 


-0.11 


0.39 


(Exp.) 


11 


-0.21 


0.00 


0.54 


RMS (TD-DFT) 


13 




0.32 


0.78 


(TBE) 


13 


0.38 


0.18 


0.48 


(Exp.) 


11 


0.25 


0.21 


0.61 


Group 2 










MSD (TD-DFT) 


36 




0.28 


0.55 


(TBE) 


36 


-0.52 


-0.24 


0.03 


(Exp.) 


12 


-0.19 


0.26 


0.59 


RMS (TD-DFT) 


36 




0.41 


0.63 


(TBE) 


36 


0.63 


0.41 


0.45 


(Exp.) 


12 


0.33 


0.38 


0.69 


Group 3 










MSD (TD-DFT) 


14 




0.42 


0.79 


(TBE) 


14 


-0.61 


-0.20 


0.17 


(Exp.) 


7 


-0.53 


-0.13 


0.27 


RMS (TD-DFT) 


14 




0.47 


0.83 


(TBE) 


14 


0.67 


0.41 


0.47 


(Exp.) 


7 


0.57 


0.39 


0.57 



TABLE III: Mean signed difference (MSD) and root-mean-square deviation (RMS) of singlet-triplet vertical excitation energies 
(in eV) calculated within TD-DFTB and TD-DFT (along the rows) with respect to TD-DFT, TBE and experiment (along the 
columns) where possible. 



TD-DFT (0.38 eV). A similar accuracy for both TD-DFT and refined TD-DFTB can be seen for the first two group 
of compounds whereas for the third group, the agreement with experimental data is somewhat better for our method 
(see |HI[ ). In addition, whereas IV still indicates some overestimation of the corrected TD-DFTB results with respect 
to TD-DFT, the MSD of our refined approach compared to experiment becomes very small, which shows that the 
corrected excitation energies are scattered around the experimental values. Specifically, in the group of aromatic 
hydrocarbons and heterocycles, excitation energies are overestimated whereas there is some underestimation for the 
carbonyl compounds. Only with regard to TBE, the refined excitation energies appear to be underestimated for every 
group of compounds. 

In contrast to the observations for singlet-triplet transitions, the traditional TD-DFTB method is, on average, in 
slightly better agreement with TD-DFT. The RMS deviations of the corrected and non-corrected methods from TD- 
DFT are 0.57 eV and 0.50 eV, respectively. On the other hand, regarding the experimental references and the TBEs, 
both approaches perform with similar accuracy. The refined formalism returns overestimated excitation energies 
according to TD-DFT results, although with respect to experiment they are again spread around the reference values. 
In this case, the overestimation for the second group of compounds is compensated with the underestimation for the 
aliphatic hydrocarbons and the nucleobases, thus leading to an almost vanishing MSD (-0.01 eV). The results for the 
latter mentioned group of compounds are the least overestimated with respect to TD-DFT, with a MSD of only 0.09 
eV. By contrast, the underestimation with respect to experiment is significant (MSD = -0.37 eV) and increases even 
more by taking the TBEs as the reference, for which the MSD amounts to -0.75 eV. For this group, the RMS error 
of TD-DFTB with respect to TBEs is remarkably high (0.91 eV and 0.96 eV for the corrected and non-corrected 
approaches, respectively). The limited agreement between the TD-DFTB excitation energies of the nucleobases and 
their TBE counterparts has been already pointed out by Trani and co-workers.^ However, it is necessary to notice 
that this failure should be rather attributed to TD-PBE itself and not to TD-DFTB as an approximation. In fact, the 
worst agreement between TD-DFT and TBE along the benchmark set is found for the carbonyl compounds and the 
nucleobases, with RMS errors of 0.96 eV and 0.89 eV, respectively. A very recent study by Foster and Wong indeed 
shows that conventional semi-local functionals fail in the description of the optical properties of nucleobases, while 
tuned range-separated functionals offer significant improvements.^ 
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count TD-DFT TD-DFTB 









new 


old 


Singlets 










MSD (TD-DFT) 


104 




0.25 


0.14 


(TBE) 


104 


-0.55 


-0.29 


-0.40 


(Exp.) 


80 


-0.28 


-0.01 


-0.12 


RMS (TD-DFT) 


104 




0.57 


0.50 


(TBE) 


104 


0.71 


0.63 


0.66 


(Exp.) 


80 


0.56 


0.48 


0.47 


Triplets 










MSD (TD-DFT) 


63 




0.29 


0.63 


(TBE) 


63 


-0.50 


-0.20 


0.14 


(Exp.) 


30 


-0.27 


0.07 


0.50 


RMS (TD-DFT) 


63 




0.40 


0.71 


(TBE) 


63 


0.59 


0.38 


0.46 


(Exp.) 


30 


0.38 


0.33 


0.63 


Total 










MSD (TD-DFT) 


167 




0.27 


0.33 


(TBE) 


167 


-0.53 


-0.26 


-0.20 


(Exp.) 


110 


-0.28 


0.01 


0.05 


RMS (TD-DFT) 


167 




0.51 


0.59 


(TBE) 


167 


0.67 


0.55 


0.59 


(Exp.) 


110 


0.52 


0.44 


0.52 



TABLE IV: Global statistics for vertical excitation energies (in eV) calculated within TD-DFTB and TD-DFT (along the rows) 
with respect to TD-DFT, TBE and experiment (along the columns) where possible. 



The overall analysis (see bottom of IV ) leads to somewhat smaller RMS errors for corrected TD-DFTB compared 
to the non-corrected formalism. It should be indicated, that despite the important improvements for the triplets 
states within the corrected formulation, the benchmark set for singlet-singlet excitations is comparatively larger, 
conceding more importance to these kind of transitions within the total balance. As a general behavior, it can be 
stated that TD-DFTB singlet excited states are shifted up in energy when the onsite correction is switched on. 
Inversely, our correction tends to shift triplet states down. This effect is particularly noticeable for n — > tt* and 
a — > 7r* excitations. The excitation energies for these transitions are either identically equal or (for few cases) very 
close to their corresponding KS energy differences at the non-corrected TD-DFTB level, resulting in degenerate or 
nearly-degenerate singlet and triplet states (see Supporting Information). The only exception for this statement are 
the triplet B\ u states of tetrazine, for which an energy displacement occurs also for traditional TD-DFTB, although 
clearly underestimated with respect to the TD-DFT results. Within the onsite correction the excitation energies 
are either shifted down for triplets or shifted up for singlets with respect to the KS energy difference, leading to 
a singlet-triplet gap in accordance with the observations at the TD-DFT level. A similar degeneracy breaking was 
shown earlier for Ng. 

Considering only n — > ir* and a — > tt* transitions, the RMS error of singlet-triplet excitation energies at the corrected 
(non-corrected) TD-DFTB level as compared to TD-DFT is 0.58 (0.82) eV. These kind of excitations are evidently 
difficult cases for the traditional formalism and an important improvement is obtained with the onsite correction, 
although the major quantitative enhancement of the latter approach is for singlet-triplet tt — > n* transitions. Indeed, 
n — >• tt* and a — > tt* excitation energies are still strongly overestimated at the corrected TD-DFTB level with respect 
to TD-DFT, with MSD of 0.49 and 0.32 eV for triplet and singlet states, respectively. However, if we compare our 
findings with the TBEs, those are by contrast, significantly underestimated (MSD = -0.21 and -0.40 eV for singlet- 
triplet and singlet-singlet transitions, respectively) and the RMS deviations for both corrected and non-corrected 
formalisms are similar. 
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TD-DFTB (old) 



10 



TD-DFTB (new) 



TD-PBE (TZP) 




2040 







700 



1,400 



2,100 



CPU time (s) 



FIG. 2: Total computational time required for calculating the low-lying excitation energies of a set of 28 molecules with 
traditional TD-DFTB (old), refined TD-DFTB (new) and full TD-DFT (PBE/TZP). 



We have introduced and implemented new extensions of the time-dependent density functional based tight binding 
(TD-DFTB) method. By taking the formalism beyond a mere Mulliken approximation, we have been able to surpass 
the wrong description of n — > ir* and a — > 7r* excitations within TD-DFTB. It was shown that particularly for the 
molecules O2, NO and N 2 these kind of transitions play an important role in the low energy optical spectra. Our 
refined formulation is in this case essential to obtain a qualitatively correct spectrum. Moreover, for the close-shell 
molecule N2 , the wrong a n* triplet-singlet degeneracy encountered by using the original approach was successfully 
overcome. We additionally identified the erroneous degeneracy of 7r — ¥ n* excitations with irreducible representations 
£~ and A for the investigated diatomic molecules as another failure of traditional TD-DFTB. The new formulation 
now delivers the correct behavior as compared to TD-DFT observations. Along with the qualitative improvement, a 
better numerical agreement with TD-DFT results can also be perceived for these systems. 

To extend the analysis of the performance of our correction, we calculated the low-lying vertical excitation energies 
for a set of 28 benchmark compounds. We found a considerable improvement in the agreement with TD-DFT and 
experiment in terms of singlet-triplet excitation energies, whereas singlet-singlet energies have the same overall quality 
as in the traditional scheme. The accuracy of the corrected formalism is similar to that of TD-DFT for both triplet 
and singlet states, whereas the computational time of the former method is roughly eighty times smaller (seeKj). 



Financial support from the Deutsche Forschungsgemeinschaft (GRK 1570 and SPP 1243) and the DAAD is gratefully 
acknowledged. 



Let {<p a } and {4>p} be a complete set of real orbitals placed at atom A and B, respectively. The orbitals are assumed 
to be orthonormal within each set, while the overlap between orbitals on different atoms is given by S a p. Let also 
f£B, kgC and Aefl unless otherwise specified. Then, the orbitals </> M and </>„ can be expanded as 
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Appendix A: Onsite correction to the total energy 




(Al) 



a<EA 



Hence, the orbital product can be expressed as 




(A2) 
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and more conveniently as 

^(r)^(r) = l -S^ (|^(r)| 2 + |^(r)| 2 ) + ± ( £ S^ Q (r)<^(r) + £ Sp^WM*) | , (A3) 
where the first term accounts for the Mulliken approach. Let us now denote with (^lv\kX) a two-electron integral with 



1 (o&p. &±v \ 



an arbitrary kernel. Using A3 the integrals (/iz^KA) can then be expanded as 



1 

4 

'l^K s^x 



{fiv\n\) = -S^ U S K \ [(fi/j,\KK) + (nn\\\) + {vv\kk) + (vu\XX)] 



( E 5 7A [(Wl«7) + (H«7)] + E ^ KM^) + (w|<5A)] 
^S KX S »- [(H«k) + (HAA)] + E V K^l««) + (fiv\XX)] 

\a£A 0£B 
^ I a^M 7#fi "t^M <5^A \ 

2 [ E E S av S~,\(pa\Krf) + E E SavStitipaltiX) 

7 EE S^S lX (f3u\K>l) + E E SpM^\SX) . (A4) 



4 



4 



Note that the expression A4 is exact as long as the AO sets {4>a}, {4>fi}i {^7} an( A {4>d} are complete. The first 
line in eq |A4| contains the leading terms, which include Coulomb-like integrals. Truncation of the expansion up to 
the first line accounts for the Mulliken approach. A next level of approximation demands the inclusion of fully-onsite 
exchange-like integrals, i.e. {fii/\fiv), with /i,ce4 and fi ^ v. At this level of theory the general four-center integrals 
can be expressed as 



{pbv\nX) « {Hv\K\)man + (/H"^)ons, ( A5 ) 



where 



(fiv\nX) m un = jS^Skx [(fifj,\KK) + (mmI^A) + (w\kk) + (vv\XX)\ , (A6) 



and 



(fiv\KX) oaB = ^ E S a vS a \(fia\(j,a)6,j. K + ^ E S aiyS aK {^a\^a)S tJi x 

a£A aeA 



• / j ~ pi*— ys\\r- - 11-' I'm ' , 

+ ^<S , ^S , AlA (AiK|MK)^Ac(l ~ <W) + j^A^reGuAl/xA^rKl - <^a) 
+ ^S^S u \(vk\vk)5 B c{^ ~ 6vk) + ^S\ l ^S l/K (uX\uX)6 BD (l - S yX ) ■ (A7) 
The contribution to the DFTB total energy originating from this correction is written as 

= \ E E AP ^ Hfc ^o]l«A) ons AP; A . (A8) 

err iivk\ 
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After substituting AT in ~KS we finally arrive at 



E on S = EE E Ap ^ M/£s>]M ap;„, (A9) 

ar A fj,v£A 

where P° u are the matrix elements of the dual density matrix defined in|8j 

Appendix B: Onsite correction to dipole matrix elements 

To derive an expression for the KS dipole matrix elements, we expand the KS orbitals into a set of localized 
atom-centered AO, ip sr7 = Y]„ C JL^V- Thus, the dipole matrix elements read 

Wwifi^> = E c ^i f i^< t ( B1 ) 

Let n £ A and v G B unless otherwise indicated. Using the orbital product expansion |A3| the AO dipole matrix 
elements, can be expressed as follows, 

(/i|r» = -V (R A + K B ) + - E s °M\i Im) + E S M*\v) , (B2) 

\a£A fieB J 



where = (/i|f \fi) and = (v\r \v) denote the positions of centers A and B, respectively. After substituting B2 
in |Bl[ we finally have 

(^if|^> = E R ^+E E m a (^)> ( B3 ) 

A A nueA 

where the definition |10| was additionally employed. 
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